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1.0  ABSTRACT 

A  previous  report  described  a  higher-order  three-dimensional  panel  method 
that  represents  the  body  about  which  flow  Is  to  be  computed  by  means  of  curved 
four-sided  surface  panels  having  linearly  varying  source  and  vorticity  distri¬ 
butions.  The  method  of  accounting  for  lift  was  incomplete*  and  the  main  purpose 
of  the  present  work  was  to  remedy  this  defect.  A  number  of  other  modifications 
to  the  method  were  also  made  to  improve  its  efficiency  and  accuracy. 


The  modifications  documented  here  are:  formulas  for  the  edge-vortex 
influences,  which  previously  had  been  neglected;  new  surface  vorticity  formulas 
that  express  its  influences  in  terms  of  source  influences;  a  modified  global 
vorticity  algorithm  to  improve  continuity  over  the  surface,  and,  an  extrapolated 
Kutta  condition. 
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3.0  PRINCIPAL  NOTATION 


area  of  the  flat  projected  panel  (Figure  1) 

derivative  of  the  equivalent  dipole  strength  along  an  N-line 

constant  governing  the  equivalent  dipole  term  quadratic  in  n 

used  as  subscripts  to  denote  quantities  associated  with  the 
first  and  second  N-lines  of  a  panel,  respectively 

arc  length  along  an  N-line  from  the  trailing  edge  to  the  n-axis 
of  panel  coordinates 

unit  vectors  along  the  axes  of  the  panel  coordinate  system 

functions  defined  by  equations  (6.3.1),  (6.3.2)  and  (6.3.4) 

total  arc  length  of  an  N-line  from  lower  surface  trailing  edge 
to  upper  surface  trailing  edge 

panel  "curvatures,"  i.e.,  second  derivatives  of  the  shape  of 
the  curved  panel  at  the  origin  of  panel  coordinates 

magnitude  of  r 

vector  from  a  point  on  the  curved  panel  to  a  point  in  space 
where  velocity  is  to  be  evaluated 

vector  from  a  point  on  the  projected  flat  panel  to  a  point  in 
space  where  velocity  is  to  be  evaluated 

surface  area  of  the  curved  panel  ( Figure  1) 

slope  of  an  N-line,  eq.  (6.2.3) 

velocity  vector 

parameters  in  the  parametric  cubic  fit 
width  of  a  panel  between  N-lines 

coordinates  of  a  point  in  space  where  velocity  is  to  be 
evaluated 

quadratic  surface  approximating  the  curved  panel,  eq.  (5.3.2) 


coordinates  of  the  first  and  second  N-lines,  respectively,  in 
panel  coordinates 

equivalent  dipole  strength.  When  used  with  subscripts  x  and 
y  it  denotes  the  corresponding  derivative  of  u  at  the  origin 
of  panel  coordinates 
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coordinates  of  a  point  on  the  curved  panel  in  panel  coordinates. 
Setting  c  =  0  gives  the  corresponding  point  on  the  projected 
flat  panel 

source  density 

vorticity  vector,  used  with  various  subscripts 
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4.0  INTRODUCTION 

Reference  1  describes  a  so-called  first-order  panel  method  for  calculating 
potential  flow  about  arbitrary  three-dimensional  lifting  bodies.  The  program 
deck  corresponding  to  this  method  has  been  distributed  widely  and  Indeed  has 
become  a  standard  design  tool  for  subsonic  aerodynamic  analysis  throughout  the 
country.  Recently  a  higher-order  version  of  the  panel  method  was  constructed 
as  described  in  reference  2.  Compared  to  the  first-order  method  of  reference  1, 
the  higher-order  method  achieves  greater  calculational  accuracy  for  a  given 
panel  number  and  thus  has  the  potential  to  reduce  cost  substantially  for  a  given 
accuracy  and  to  calculate  flow  about  more  complicated  configurations.  The 
initial  nonlifting  version  of  the  higher-order  method,  which  was  documented  in 
reference  3,  has  already  been  used  in  design  applications,  for  example  refer¬ 
ence  4.  However,  complete  formulas  and  logical  procedures  for  the  method  were 
first  presented  in  reference  2,  which  also  included  provisions  for  lifting 
effects.  The  lifting  method  described  in  reference  2  to  a  large  extent  has 
the  nature  of  a  pilot  method  for  proving  the  approach,  rather  than  a  final 
general  procedure. 

The  work  described  in  this  report  consists  of  modifications  to  the  method 
of  reference  2  that  greatly  increase  its  accuracy  and  numerical  efficiency. 

To  make  the  present  document  self-contained  would  require  a  complete  description 
of  the  higher-order  panel  method  and  thus  a  duplication  of  large  portions  of 
reference  2.  Since  reference  2  is  generally  available,  the  consequent  large 
increase  in  the  length  of  the  present  report  seems  neither  necessary  nor 
desirable.  Accordingly,  reference  2  is  relied  upon  to  provide  the  general  ideas 
and  philosophy  of  the  higher-order  method,  as  well  as  detailed  formulas  and 
logic.  Indeed  references  will  frequently  be  made  to  sections  and  even  indiv¬ 
idual  equations  of  reference  2.  A  few  features  of  the  higher-order  method 
that  bear  directly  on  the  modifications  described  in  this  report  are  mentioned 
below.  But  the  intention  is  one  of  emphasis  not  of  completeness. 

A  panel  method  discretizes  the  body  about  which  flow  is  to  be  computed, 
representing  it  by  a  large  number  of  small  four-sided  surface  panels,  on 
each  of  which  are  distributions  of  source  and  dipole  or  vortlcity.  A  control 
point  is  selected  on  each  panel  where  the  normal -velocity  boundary  condition 
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is  enforced  and  where  flow  quantities  are  eventually  calculated.  In  lifting 
cases  a  Kutta  condition  is  applied  to  insure  smooth  flow  off  wing  trailing 
edges.  In  essence  the  source  strengths  on  the  panels  are  adjusted  to  satisfy 
the  normal-velocity  boundary  conditions,  and  thus  each  panel  has  an  independent 
value  of  source  strength.  The  variation  of  dipole  or  vorticity  strength  over 
the  surface  Is  fitted  by  certain  global  algorithms  to  obtain  expressions  con¬ 
taining  a  number  of  adjustable  parameters  equal  to  the  number  of  locations 
where  the  Kutta  condition  is  applied. 

The  key  calculational  unit  in  a  panel  method  is  the  set  of  formulas 
giving  the  velocities  induced  at  a  point  in  space  by  the  source  and  dipole  or 
vorticity  distributions  on  a  panel.  In  the  first-order  method  of  reference  1 
the  panels  are  planar  and  the  required  induced  velocities  due  to  constant 
source  and  quadratic  dipole  distributions  may  be  obtained  exactly  by  means  of 
analytic  integration  over  the  panel.  When  the  point  in  question  is  far  from 
the  panel,  approximate  "far- field"  formulas  are  used  to  reduce  computing  time. 
In  the  higher-order  method  of  reference  2,  the  panels  are  conceptually  curved. 
Since  analytic  integration  over  a  curved  panel  is  not  possible,  the  induced 
velocities  due  to  linearly  varying  source  and  vorticity  distributions  on  the 
curved  panel  are  expressed  as  expansions  about  the  effects  of  the  flat  panel 
that  is  the  projection  of  the  curved  panel  in  the  tangent  plane.  Integration 
over  the  "projected  flat  panel"  of  the  various  terms  in  the  expansion  can  be 
performed  analytically.  The  validity  and  consistency  of  the  expansions  are 
discussed  in  reference  2.  Approximate  expressions  for  the  terms  of  the 
expansion  are  used  at  distant  points.  The  errors  due  to  using  these  far-field 
expressions  are  independent  of  and  can  be  made  small  with  respect  to  the 
truncation  error  of  the  expansion  about  the  projected  flat  panel. 

The  analytic  expressions  in  reference  2  for  the  effects  of  panel  vorticity 
are  extremely  cumbersome.  Moreover,  they  do  not  lend  themselves  to  efficient 
far-field  approximations.  The  result  is  long  computing  times  for  the  method 
as  presented  In  reference  2.  This  situation  has  been  remedied  by  relating 
the  vorticity- Induced  velocities  to  the  source- induced  velocities  and  thus 
obtaining  the  former  In  essentially  no  additional  time.  The  far-field  problem 
no  longer  arises  because  the  source  far-field  approximation  is  all  that  is 
required.  This  modification  is  discussed  in  section  5.0. 
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As  discussed  in  reference  2,  vorticity  has  been  employed  on  a  panel  rather 
than  a  dipole  distribution  because  it  leads  to  a  simpler  expansion  about  the 
projected  flat  panel.  Thus  the  effects  of  line  vortices  around  the  perimeter 
of  the  panel  are  neglected.  Over  the  portions  of  the  surface  away  from  any 
physical  edges  the  effect  of  the  line  vortex  along  each  edge  of  a  panel  is 
cancelled,  at  least  to  some  order,  by  that  of  the  line  vortex  along  the  edge  of 
the  adjacent  panel.  The  result  is  a  relatively  weak  line  vortex.  At  physical 
edges  of  the  body,  e.g. ,  a  wing  tip,  there  can  be  no  such  cancellation,  and 
indeed  a  strong  edge  or  tip  vortex  is  known  to  be  present.  It  turns  out  that 
it  is  the  trailing  edge  vortices  that  are  important.  Section  6.0  presents  form¬ 
ulas  that  account  for  the  effects  of  such  edge  vortices. 

So-called  global  vorticity  algorithms  are  used  to  relate  the  vorticity 
strengths  on  the  panels  and  thus  reduce  the  variation  of  vorticity  over  the 
surface  to  analytic  expressions  depending  on  a  number  of  parameters  equal  to 
the  number  of  locations  at  which  the  Kutta  condition  is  prescribed.  Section  7.0 
describes  two  improvements  to  the  global  vorticity  algorithms  of  reference  2. 

One  is  concerned  with  minimizing  spanwise  vortices  at  interior  panel  edges 
(previous  paragraph)  and  the  other  with  an  improvement  in  the  spanwise  vorticity 
fit  that  eliminates  extraneous  wake  vorticity. 

The  Kutta  condition  used  in  all  versions  of  the  panel  method  is  an  equal - 
pressure  condition  applied  at  the  upper  and  lower  trailing  edge  of  a  wing.  In 
references  1  and  2  the  pressures  set  equal  are  those  at  the  uDDer  and  lower 
control  points  nearest  the  trailing  edge  for  each  spanwise  location.  Thus  the 
equal -pressure  condition  is  applied  at  a  distance  of  half  a  panel  length  from 
the  trailing  edge.  This  is  an  acceptable  approximation  for  a  first-order 
method,  but  it  penalizes  a  higher-order  method  unnecessarily.  Section  8.0 
describes  an  Improved  Kutta  condition  that  extrapolates  upper  and  lower  surface 
pressures  to  the  trailing  edge  before  setting  them  equal. 

finally,  the  geometric  procedure  described  In  reference  2  for  calculating 
panel  curvatures,  control  points,  and  normal  vectors  has  been  found  to  be 
unsatisfactory  (in  some  cases).  Accordingly,  it  has  been  replaced  by  a  more 
elaborate  procedure.  Although  this  modification  is  outside  the  scope  of  the 
work  reported  here.  It  seemed  that  for  completeness  a  description  of  it  should 
be  included.  Thus  it  is  presented  in  Appendix  A. 
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5.0  CALCULATION  OF  VORTICITY  INDUCED  VELOCITIES  IN  TERMS  OF 
SOURCE  INDUCED  VELOCITIES 


5.1  Background 

As  has  been  mentioned,  the  expressions  for  induced  velocity  due  to  any 
polynomial  distribution  of  singularity  can  be  analytically  integrated  over  a 
flat  panel.  Moreover,  the  only  trancedental  functions  that  occur  are  the 
logarithms  and  inverse  tangents  that  are  already  present  in  the  first-order 
constant-source  formulas.  Thus,  the  expressions  for  higher-order  effects,  which 
are  obtained  by  integration  over  the  projected  flat  panel,  can,  at  least 
potentially, be  evaluated  in  computing  times  only  modestly  larger  than  those  for 
the  first-order  formulas.  This  was  certainly  true  for  the  two-dimensional 
higher-order  method  of  reference  5  and  to  a  lesser  extent  for  the  nonlifting 
version  of  the  present  method  (reference  3).  However,  there  are  limits  to 
the  insensitivity  of  computing  time.  The  formulas  for  the  vorticity  effects 
that  are  presented  in  section  6.3  of  reference  2  are  so  elaborate  that  they 
add  substantially  to  the  required  computing  time.  Most  of  the  complications 
arise  from  the  curvature- dependent  terms.  Moreover,  these  formulas  are  not 
suitable  for  efficient  far- field  approximation.  Thus  the  90%  or  so  of  the 
velocity  influences  that  are  calculated  by  far-field  formulas  require  consid¬ 
erably  more  computing  time  for  the  vorticity  effects  than  for  the  source  effects, 
the  latter  of  which  do  have  efficient  far-field  approximations. 


5.2  General  Theory 

The  calculation  of  the  vorticity  influences  can  be  made  much  more  efficient 
by  expressing  them  in  terms  of  the  corresponding  source  influences,  which  of 
course  must  be  calculated  in  any  event.  The  use  of  this  procedure  was  put 
forward  in  reference  6.  The  portion  of  the  theory  that  Is  needed  for  the  present 
purpose  is  quite  easy  to  state. 


Suppose  there  is  a  variable  source  density  a  on  a  portion  of  a  plane  or 
curved  surface  S.  The  velocity  due  to  this  at  a  point  (x,y,z)  is 


$  (source)  •//VdS 


(5.2.1) 
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where  (xq,  yq,  zq)  is  a  point  on  s  and  where 

r  *  (x  -  xq)T  +  (y  -  yq)j  +  (z  -  zq)£  (5.2.2) 

As  usual  r  is  the  magnitude  of  r  (see  Figure  1).  If  there  is  a  vorticity 
distribution  of  S  of  strength 

u  =  at  (5.2.3) 

the  Biot-Savart  law  gives  the  resulting  induced  velocity  as 

$  (vorticity)  =  fj  mdS  (5.2.4) 

S  r 

Then  if  t  is  a  constant  vector  and  if  i»  has  the  same  spatial  variation 
as  a,  the  velocity  due  to  the  vorticity  distribution  may  be  expressed  in 
terms  of  the  velocity  due  to  the  source  distribution  as 

$  (vorticity)  stxl  (source)  (5.2.5) 

since  u  can  be  resolved  into  components,  each  of  which  has  a  constant  direc¬ 
tion,  the  restriction  to  a  constant  ?  is  not  serious.  Although  the  above 
results  apply  to  a  curved  surface  S,  it  is  far  simpler  to  apply  to  a  flat 
surface.  In  the  present  context  the  above  is  applied  to  the  flat  projected 
panel . 


5.3  Calculation  of  the  Velocity  Induced  by  a  Panel 

Figure  1  illustrates  the  projection  of  a  curved  panel  S  on  the  surface 
to  a  flat  panel  A  in  the  tangent  plane.  In  particular.  Figure  1  illustrates 
r  as  given  in  (5.2.2)  above  and  the  vector  r^  from  a  point  of  the  projected 
flat  panel  to  the  point  (x,y,z) 

=  (x  -  c)te  +  (y  -  n)Je  +  zl<e  (5.3.1) 

where  t  ,  j  ,  £  are  unit  vectors  along  the  axes  of  the  panel  coordinate 
system.  The  vertical  distance  c  between  the  curved  panel  and  its  projection 
is  approximated  by  its  leading  term  w^ch  represents  a  surface  of 

second  degree 
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S2  =  P?2  +  2QCn  +  Rn2  (5.3.2) 

where  P,  Q,  R  are  the  second  derivatives  of  c  at  the  origin  of  panel 
coordinates  —  the  so-called  surface  curvatures. 


The  aim  is  to  parallel  the  development  of  section  5.3  of  reference  2  and 
express  the  results  in  terms  of  source  effects.  From  equation  (53)  of  refer¬ 
ence  2  it  is  seen  that  a  two-term  expansion  of  the  vector  vorticity  distribu¬ 
tion  is 


(5.3.3) 

(5.3.4) 


“1  =  2(uxyC  +  V^e  “  2(wxxc  +  V^e  + 

+  2[-(Qc  +  Rn)wx  +  (Pc  +  Qn)uy]te  (5.3.5) 

is  first  order.  The  constants  ux»vxx.  etc.  are  the  derivatives  of  the 
equivalent  dipole  distribution  as  given  by  equation  (50)  of  reference  2. 

From  (5.3.1)  and  (5.2.2)  expressed  in  panel  coordinates: 

r  =  rf  -  (5.3.6) 

Thus  a  two-term  expansion  of  the  velocity  at  (x,y,z)  due  to  the  vorticity 
on  the  panel  is 


By  taking  the  gradient  of  the  and  terms  of  the  source  expan¬ 

sion  on  page  20  of  reference  2,  it  can  be  seen  that  the  integral  multiplying 
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above  is  just  the  sum  of  superscript  0  and  c  source  terms  for  unit 
source  density.  Specifically  this  combination  is  the  velocity 

=  $(°)  +  +  2q$(Q)  +  tfWj  (5.3.9) 

where  the  quantities  on  the  right  of  (5.3.9)  are  as  defined  in  section  6.2  of 
reference  2,  and  the  combined  velocity  V'  is  explicitly  calculated  by  the 
existing  code. 

To  analyze  the  remaining  term  of  (5.3.8),  collect  terms  in  equation  (5.3.5) 
to  obtain 


where  the  vectors 


=  2(cqx  +  nQy) 


a  =  u  1  —  u  l 

Mx  Mxy  e  MxxJe 

*>  ■+ 
q  =  u  i  —  u  i 

y  yy  e  xyJe 


+  -  0Mx)te 

+  (Quy  -  Ryx)lce 


(5.3.10) 


(5.3.11) 


are  constants  in  the  integration.  The  integrals  that  result  from  using 
(5.3.10)  in  the  second  term  of  (5.3.8)  are  the  velocities  due  to  linearly 
varying  source  densities  in  the  £  and  n  directions  having  unit  slope,  i.e. 
V(1x)  and  V(ly)  from  section  6.2  of  reference  2. 

Thus  the  velocity  (5.3.7)  due  to  vorticity  on  the  panel  may  be  expressed 
in  terms  of  source  velocities  as  follows 

Vu  =  Z0  x  V'  +  2[qx  X  V(1x)  +  qy  x  $(ly)]  (5.3.12) 

It  is  interesting  to  note  that  (5.3.12)  can  be  evaluated  directly  in  reference 
coordinates  after  the  relevant  source  velocities  have  been  calculated  and  put 
into  this  system.  As  regards  the  velocities  due  to  the  vorticity,  this  not 
only  means  that  no  transformations  between  panel  and  reference  coordinates 
are  required,  but  it  also  means  that  the  question  of  far-field  calculation 
need  never  arise.  If  the  source  velocities  have  been  computed  by  far-field 
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formulas,  they  simply  are  used  in  (5.3.12),  so  that  in  effect  the  vorticity 
calculation  uses  the  source  far-field  procedure.  The  present  code  takes 
advantage  of  the  second  of  these  facts,  the  use  of  far-field  source  formulas, 
but  performs  the  calculation  in  panel  coordinates. 

Formulas  (5.3.4),  (5.3.11)  and  (5.3.12)  replace  the  elaborate  formulas 
of  section  6.3  of  reference  2. 

5.4  Assembly  of  the  Vorticity  Formulas 

The  formulas  of  section  5.3  give  the  difficult  portion  of  the  vorticity 
calculation.  However,  the  task  remains  to  collect  the  terms  of  the  vorticity- 
induced  velocity  into  velocities  associated  with  the  first  and  second  N-lines 
of  the  panel  in  a  manner  similar  to  that  of  section  9.0  of  reference  2. 
Basically  this  is  a  matter  of  performing  the  indicated  vector  cross  products 
of  section  5.3  above  and  using  formulas  (141)  of  reference  2  for  the  y 
derivatives.  This  last  is  most  easily  done  computationally  by  means  of  the 
following  logical  table 


u-derivative 


First  N-line 

h  w 
/-c(nl  +  n3} 

d 

2w 


Second  N-line 


"i T+  c(nl  +  n3} 

-d  (5.4.1) 

1_ 

2w 


Pyy 


c 


-c 


where  all  quantities  have  the  same  meaning  as  in  reference  2  section  9.2.  The 
foregoing  are  on-body  formulas.  For  wake  panels  set 


*  u 


XX 


=  0 


hp  =  Lp  (total ) , 


h$  =  Ls  (total) 


(5.4.2) 
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In  performing  the  indicated  cross  products  it  will  be  recalled  that  the 

components  of  Tg*  Je»  (k  *  n)  in  reference  coordinates  are  the  entries 

of  the  panel's  transformation  matrix  a„„,  i.e. 

mn 

Te  ’  allT  *  a12*  +  “13* 

Je  *  a21t  +  a22J  +  a23?c  (5.4.3) 

*e  =  a31T  +  a32J  +  a33* 

If  the  calculation  is  performed  in  panel  coordinates,  e.g.  the  near  field,  these 
are  simply  the  unit  vectors  along  the  coordinate  axes,  and,  for  example  the 
coefficient  of  T  is  simply  V  in  the  panel  coordinates. 

6  A 

The  formulas  of  section  5.3  above  and  the  present  section  permit  deter¬ 
mination  of  the  F  and  S  vorticity  influences  in  the  form  of  equation  (143) 
of  reference  2  in  a  straightforward  way. 
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6.0  THE  LINE  VORTEX  ALONG  A  STREAMWISE  EDGE  OF  A  PANEL 

6.1  Background 

As  mentioned  in  the  Introduction  and  discussed  in  section  5.3  of  reference 
2,  the  use  of  a  vorticity  distribution  on  a  panel  rather  than  a  dipole  distri¬ 
bution  means  that  the  effects  of  a  line  vortex  around  the  perimeter  of  the  panel 
are  neglected.  Away  from  a  physical  edge  of  the  body  this  neglect  might  be 
justified  because  line  vortices  on  the  edges  of  adjacent  panels  tend  to  cancel. 
It  turns  out  this  is  true  for  s panwise  panel  edges  but  not  for  streamwise  edges. 
That  is,  the  edge-vortex  contribution  to  the  bound  vorticity  has  been  found  to 
be  relatively  unimportant  but  the  trailing  edge  vorticity  along  the  N-lines  must 
be  accounted  for.  This  matter  is  discussed  more  in  section  7.0  and  examples  are 
given  in  section  8.0.  At  a  wing  tip  or  other  physical  edge  of  a  lifting  body 
there  is  a  strong  edge  vortex  whose  neglect  renders  the  solution  significantly 
nonpotential.  The  absence  of  machinery  for  calculating  the  effects  of  such  an 
edge  vortex  restricts  the  method  of  reference  2  to  very  specialized  lifting 
configurations.  The  analysis  of  this  section  removes  this  restriction.  The 
same  formulas  are  used  for  interior  edge  vortices  along  the  N-lines. 

There  are  several  ways  of  accounting  for  the  effect  of  the  edge  vortex,  all 
of  which  are  theoretically  equivalent  to  some  order  of  accuracy.  The  approach 
used  here  is  the  analogy  of  that  used  throughout  the  higher-order  development. 

A  vortex  lying  along  the  edge  of  the  curved  panel  is  projected  into  the  tangent 
plane. 

6.2  Derivation  of  the  Influence  of  an  Edge  Vortex 

The  equation  of  the  curved  panel  is  (5.3.2).  For  definiteness  consider 
the  case  when  the  edge  in  question  lies  in  the  plane  n  ■  n-j ,  i.e.  the  first 
N-line.  The  modifications  for  the  case  of  the  second  N-line  are  obvious. 

Thus  the  curve  c  along  which  the  vortex  lies  is 

C  =  c(c)  -  PC2  +  2Q5i>j  +  Rn2  (6.2.1) 

The  unit  vector  along  this  curve  is 

t  -  ■y-1--  [t  +  Tj  1  (6.2.2) 

VTT7 
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where 

T  =  Z{H  +  Qn-j ) 

The  velocity  due  to  the  vortex  is 


(6.2.3) 


?.  -  Wh  (6.2.4) 

c  r 

where  u  is  the  edge  value  of  the  equivalent  dipole  strength.  Arc  length 
along  the  curve  is  related  to  distance  in  the  tangent  plane  by 

ds  =^1  +  T2  d£  (6.2.5) 

Thus  with  r  expressed  in  panel  coordinates  (Figure  1) 


(t  x  r)ds  =  {[  0 

-  T(y  —  n-| )  ] 

Te 

[  z 

+  T(x  -  0  +  t  1 

Je 

(6.2.6) 

[(y  -  n-,) 

+  0  ] 

where  the  terms  in  the  first  column  of  (6.2.6)  are  first  order,  and  those  in 
the  second  column  are  second  order.  This  expression  is  exact  except  for  the 
approximation  c  =  C2. 

3 

As  shown  in  reference  7  a  three  term  expansion  of  1/r  is 


=  [1  +  3c-j  +  3(c^  +  c^)] 

r  rf 


(6.2.7) 


where 


C1  = 


r  -  3  2  _  1  c2 

c2  '  2  C1  I  “? 


(6.2.8) 


Along  the  N-line  the  equivalent  dipole  strength  varies  linearly 

v  s  B(h  +  e) 


(6.2.9) 


where  h  is  the  total  arc  length  along  the  N-line  up  to  the  n-axls  of  panel 
coordinates  (see  section  9.2  of  reference  2),  and  B  is  the  unknown  value  of 
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vorticity  that  is  determined  from  the  Kutta  condition.  The  fundamental  flow 
is  obtained  by  setting  B  equal  to  unity  1*n  (6.2.9).  Multiplying  the  above 
expansions  gives  the  components  of  the  vortex  velocity  as  follows 


[T(y  -  ^)h]  -  [T(y  -  ^JOc^h  +  C)] 


d? 


52 

Vry  ■  /  ^  j-zh  +  (-z(3Clh  +  c)  +  h(T(x  -  c)  +  c2)]  (6.2.10) 

5-,  rf 

+  [-z{3(c2  +  c2)h  +  3c15>+  (3c.,h  +  c)(T(x  -  O  +  c2)l{dC 


6.3  Formulas  for  the  Edge- Vortex  Influence 
The  integrals  in  (6.2.10)  have  the  form 

Jmn=  1*$* 

<1  f 

Once  the  JQn  and  Jln  have  been  calculated  the  others  are  calculated  from 
the  recursion  formulas 

Jmn  =  J(m-2)(n-2)  +  2xJ(m-l)n  “  p2j(m-2)n  (6.3.2) 

where 

p2  =  x2  +  (y  -  n-| )2  +  z2  (6.3.3) 

The  required  JQn  and  J^n  are 

=  log  1  2  _  ■  j1 2  =  -L^12^  (eq.  (7.7.3)  reference  1) 

ui  +  r2  o12 

J11  =  r2  “  rl  +  xJ01 
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J 

J 

J 

J 

J 

J 


,  1 

03 

'7 

_  1 

13 

rl 

_  1 

C2  “  x 

X 

1 

Uf 

1 

r2 

rl 

05 '17 


—  +  xJn0 
r2  03 

r52  -  X  C1  -  X 


1 


15  =  “  I  (“J”"?)  +  xJ05 


r2  rl 


07  =i? 


C2  -  x  ^  -  x 


1 


17  =  J  (7“7}  +  X°07 
rl  r2 


+  2J 


03 


+  4J, 


05 


(6.3.4) 


where 

q2  =  (y  —  n^2  +  z2  (6.3.5) 

and  where  r1  and  r2  are,  respectively,  distances  of  the  point  (x,y,z) 
from  the  ends  of  the  interval,  i.e., 

rk =  (x  _  ck)Z  +  (y  -*i>2  +  z2  (6.3.6) 

which  is  the  same  definition  used  in  references  1  and  2. 

In  terms  of  certain  auxiliary  functions  Fn  the  velocity  components  of 
(6.2.10)  are 

V-fy-^KhF,  +  f5j 

Vpy  =  -zhJQ^  +  [-zF2  +  hFj]  —  z[3hF^  +  Fg]  +  F^  (6.3.7) 

vrz  =  (y  _nl)lhJ03  +  F2  +  3hF4  +  F6] 

The  auxiliary  functions  for  on-body  panels  are: 
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=  2PJ13  +  2C^1J03 

?2  ~  3zh[PJ25  +  2Qn-|Ji5  +  Rn]JQ5^  +  *J13* 

F3  =  2PxJ-j3  +  2QtvjxJq3  -  PJ23  +  Rn-|J03 
C  -  i  ,2n  _  1 0 

F4  "  2  2  q7  2  q5 

Qj  =  ip2°4j  +  4PQrllJ3j  +  *2PR  +  4Q2)nlJ2j  +  4QRnlJlj  +  (6,3'8) 

F5  =  6zh[P2J33  +  3PQn-|J25  +  PRo2^  '  >  T,f 5  +  qPniq05^ 

+  {  2  Pu  2  ^  1  ■  :  j  ^  1 3  ) 

Ffi  =  {3z[PJ35  +  2Qn-|  J25  +  Rn^153> 

F7  =  3zh  l-P2J45  +  (2P2x  -  2PQn1)035  +  ePQxn^^ 

+(4Q2xn2  +  2PRxn2  +  2QRn^)J15  +  (2QRxn3  +  R2nj)J05] 

+  {-PJ33  +  2PxJ23  +  {2QxDt  +  Rn2)Ju> 

The  formulas  for  wake  panels  are  obtained  from  (6.3.8)  by  deleting  all  terms 
in  {}  and  replacing  h  by  L  (total)  (section  9.2.2,  reference  2). 

For  the  semi -infinite  last  wake  additional  changes  are  made  to  the 
formulas  (6.3.4)  for  the  Jmn  corresponding  to 

^2  -*•  »  r2  -*■  ®  ^2^r2  ^  (6.3.9) 

Furthermore  P  and  Q  are  set  equal  to  zero. 

6.4  Far-Field  Formulas 

Originally  It  had  been  thought  that  edge  vortices  would  be  required  only 
along  physical  edges  of  a  lifting  body,  such  as  wing  tips.  In  such  a  case  no 
approximate  far-field  formulas  would  be  needed,  because  the  total  number  of 
edge-vortex  influences  would  be  a  small  fraction  of  the  panel  source  Influences 
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However,  numerical  experimentation  showed  that  for  good  accuracy  edge  vortices 
are  required  along  all  streamwise  panel  edges,  i.e.  along  all  N-lines  of  a 
lifting  body.  Thus  the  total  number  of  edge- vortex  influences  is  substantially 
twice  the  number  of  source  influences,  and  simple  far-field  formulas  are  needed 
to  reduce  the  computing  time.  These  are  used  along  interior  edges.  The  nota¬ 
tion  is  the  same  as  in  section  6.3,  which  refers  to  the  first  edge  of  a  panel. 


Compute 


r  F  = 


Ci  +  Co 


2  2 
+  [y  -  nn]  +  Z 


If  (C2  “  ?! )2/r|  <  0.001,  use 

Vrx  =  -(y-ni)V 

Cl  +  Co 

Vry=  [‘Z  +  To  ^x— +  ?oir 


vrz  =  (y  -  n-,)I 


where 


I  -  (h  + 

rF 

Ci  +  Co 

To  “  2tp  -5-2— S+  Qnil 


(6.4.1) 


(6.4.2) 


(6.4.3) 


Ci  +  Co  2 

-  PC-1-*-1) 


+  2Qn^  ( 


Ci  +  c2 


T-±)  +  RnJ 


For  the  second  N-line  the  obvious  quantities  are  replaced  by  the  corresponding 
ones.  The  above  equations  replace  the  elaborate  formulas  of  the  previous 
section. 
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7.0  MODIFICATION  OF  THE  GLOBAL  VORTICITY  ALGORITHM 
7.1  Considerations  of  Continuity  and  Accuracy 

In  general  the  singularities  are  discontinuous  from  one  panel  to  the  next 
over  the  surface  of  the  body,  i.e.,  across  the  interior  panel  edges.  In  a 
consistent  method  these  discontinuities  are  higher  order  in  some  sense,  but 
it  seems  intuitively  desirable  to  reduce  or  eliminate  them.  However,  the 
very  useful  method  of  reference  1  did  not  appear  to  suffer  very  much  from  their 
presence.  The  discontinuity  in  source  strength  is  an  unavoidable  consequence 
of  the  basic  method  of  solution  as  described  in  reference  2,  section  8.0.  In 
any  event  experience  indicates  that  inaccuracies  due  to  source  discontinuities 
are  relatively  unimportant.  It  is  the  problem  of  dipole  discontinuity  that  is 
addressed  here. 

If  the  lifting  effects  are  accounted  for  by  means  of  dipole  distributions 
on  the  panels,  e.g. ,  reference  1,  then  a  discontinuity  of  dipole  strength 
yields  a  concentrated  line  vortex  along  the  edge  of  a  panel  whose  effects 
presumably  are  almost  cancelled  by  the  corresponding  edge  vortex  of  the 
adjacent  panel.  Thus  the  possible  inaccuracies  are  due  to  presence  of  weak 
line-vortex  singularities  lying  in  the  surface  where  clearly  no  such  singularity 
exists.  In  the  present  method,  however,  the  lifting  effects  are  accounted  for 
by  means  of  vorticity  distributions  on  the  panels  which  are  derived  from  an 
equivalent  dipole  distribution.  If  this  last  is  discontinuous  across  interior 
panel  edges,  there  will  be  no  line  vortex.  Instead  the  computed  velocity  field 
will  be  slightly  nonpotential.  It  is  not  obvious  which  of  the  above  two  types 
of  error  is  larger  or  if  either  is  particularly  serious.  Nevertheless,  it 
seemed  useful  to  try  to  minimize  the  effect  of  this  source  of  error. 

There  are  two  distinct  phenomena  that  cause  the  equivalent  dipole  distri¬ 
bution  to  be  discontinuous.  The  first  is  geometric.  In  general,  the  edges 
of  adjacent  panels  are  not  exactly  coincident,  and  thus  cancellation  of  edge 
singularities  cannot  be  exact.  This  situation  is  unavoidable  In  the  present 
method  and  will  not  be  considered  further.  Attention  has  been  directed  towards 
the  second  cause  of  discontinuity,  namely  nonagreement  of  the  equivalent  dipole 
distributions  of  adjacent  panels  along  their  common  edge.  There  are  two  types 
of  interior  edges  -  spanwise  and  streamwise,  the  latter  of  which  lie  along  the 


20 


NADC-79277-60 


N-lines.  Section  7.2  addresses  the  question  of  obtaining  continuity  of  the 
equivalent  dipole  distributions  across  spanwise  edges  by  suitably  modifying 
the  spanwise  dipole  variation  and  thus  the  chordwise  vorticity  variation. 
Exact  continuity  can  be  obtained  In  the  absence  of  geometric  discontinuity. 
Section  7.3  discusses  a  modification  of  the  spanwise  vorticity  variation 
that  apparently  produces  an  improvement  but  not  exact  continuity. 


7.2  Modification  of  the  Chordwise  Vorticity  on  a  Panel 


The  basic  analysis  has  been  carried  out  in  reference  1  for  the  flat  panel 
case,  and  it  applies  to  the  present  method  as  well  under  the  assumption  of 
coincident  panel  edges.  As  stated  in  section  7.3.4  of  reference  1,  the  dipole 
strength  is  in  general  discontinuous  between  two  adjacent  panels  of  the  same 
lifting  strip  unless  a  term  quadratic  in  spanwise  location  is  added  to  each 
panel  distribution.  The  form  of  this  term  is 


CAB  (n  —  n-|  )(n  —  n^) 


(7.2.1) 


where  n  is  the  spanwise  panel  coordinate,  n-j  and  n3  are  the  locations 
of  the  parallel  sides,  AB  is  the  change  in  vorticity  strength  across  the 
span  of  the  panel  (zero  if  a  piecewise  constant  B  variation  is  used)  and 
c  is  a  constant  to  be  adjusted  for  continuity.  As  shown  in  section  7.3.4 
of  reference  1,  continuity  between  panels  i  and  i  +  1  of  a  strip  is  obtained 


if 


(7.2.2) 


where  w  Is  panel  width  (usually  the  same  for  all  panels  of  a  strip),  m32 
is  the  slope  of  the  upper  panel  edge  and  m^  the  slope  of  the  lower  panel 
edge.  Equation  (7.2.2)  is  solved  for  successive  values  of  c^  beginning 
with 

c(1)  =  0  (7.2.3) 

and  proceeding  over  all  on-body  panels  of  the  strip.  The  choice  (7.2.3)  is 
arbitrary  and  expresses  the  fact  that  equations  (7.2.2)  have  a  nonunique 
solution,  once  the  values  of  c  for  on-body  panels  have  been  calculated. 

The  wake  values  are  given  by  the  procedure  of  section  7.9  of  reference  1. 
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7.3  Modification  of  the  Spanwise  Fitting  Procedure 

In  the  standard  option  of  the  present  method  the  variation  of  the  equivalent 
dipole  strength  along  an  N-line  is  Bs  where  B  is  a  constant  and  s  is  arc 
length  measured  around  the  contour  from  the  trailing  edge,  e.g.,  equation  (7.3.1) 
of  reference  1  or  (138)  of  reference  2.  Thus  B  denotes  vorticity  strength 
normal  to  an  N-line,  i.e.,  bound  vorticity.  The  total  circulation  around  the 
section  defined  by  the  N-line  is  BL  (total),  where  L  (total)  is  the  total 
arc  length  around  the  section  from  lower  trailing  edge  to  upper  trailing  edge. 

The  spanwise  fitting  option  described  in  section  7.11  of  reference 
uses  a  piecewise  constant  or  a  piecewise  linear  fit  to  the  values  of  B. 

Only  the  linear  fit  is  permitted  in  the  higher-order  method  of  reference  2. 
However,  it  is  clear  from  equation  (7.3.1)  of  reference  1  that  the  equivalent 
dipole  distribution  in  the  wake  depends  on  values  of  BL  (total)  rather  than 
B.  This  is  clear  also  on  physical  grounds  since  the  strength  of  the  trailing 
wake  vorticity  equals  the  spanwise  derivative  of  the  bound  circulation,  which 
is  BL  (total).  Accordingly  a  fit  in  terms  of  B  can  lead  to  extraneous 
wake  vorticity  that  does  not  necessarily  vanish  in  the  limit  of  small  strip 
width.  Since  the  method  of  reference  1  has  given  many  good  results,  clearly 
the  above  effect  is  not  important  in  many  cases.  Nevertheless  it  seemed  useful 
to  eliminate  this  possible  source  of  inaccuracy,  and  the  spanwise  fits  have 
been  modified  to  accomodate  BL  (total)  rather  than  B.  This  change  also 
has  the  effect  of  reducing  the  discontinuity  of  the  equivalent  dipole  distri¬ 
bution  across  N- lines  and  thus  the  magnitude  of  the  required  edge  vorti cities. 
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8.0  THE  USE  OF  AN  EXTRAPOLATED  KUTTA  CONDITION 

The  present  method  and  its  first-order  predecessor  both  use  an  equal -pressure 
Kutta  condition,  in  which  the  squares  of  the  velocities  on  the  upper  and  lower 
surfaces  of  the  trailing  edge  are  required  to  be  equal.  In  the  first-order  method 
of  reference  1  the  velocities  used  are  those  at  the  control  points  nearest  the 
trailing  edge  on  the  upper  and  lower  surfaces  of  the  wing.  Thus  in  effect  the 
Kutta  condition  is  applied  a  finite  distance,  i.e.,  half  a  panel  length,  forward 
of  the  trailing  edge.  This  approximation  is  sufficient  for  conventional  airfoils 
where  the  difference  between  the  upper-  and  lower-surface  pressures  a  short  distance 
ahead  of  the  trailing  edge  is  small.  However,  for  supercritical  airfoils  pressure 
gradients  near  the  trailing  edge  are  large  and  even  a  short  distance  forward  the 
difference  between  upper-  and  lower-surface  pressures  is  substantial.  Requiring 
these  pressures  to  be  equal  at  a  location  where  they  actually  are  quite  different 
results  in  a  considerable  loss  of  lift.  This  situation  was  obscured  for  the  first- 
order  method  by  other  inaccuracies  that  occur  in  cases  of  supercritical  airfoils. 
However,  the  higher-order  method  handles  such  airfoils  much  more  accurately  and 
the  use  of  the  older  form  of  the  Kutta  condition  penalizes  such  a  method 
unnecessarily. 

The  numerical  implementation  is  quite  simple.  Using  velocity  components  at 
the  two  control  points  nearest  the  trailing  edge  on  the  upper  and  lower  surfaces, 
linear  extrapolation  gives  values  of  these  components  "at"  the  trailing  edge.  The 
sum  of  the  squares  of  the  upper-  and  lower-surface  components  are  then  equated  in 
the  usual  way. 


NADC-79277-60 


9.0  CALCULATED  EXAMPLES 

The  formulas  presented  In  this  report  have  been  thoroughly  verified  to 
establish  that  they  are  coded  as  they  are  written.  The  larger  questions  of  the 
correctness  and  effectiveness  of  the  approach  and  the  appropriateness  of  the 
various  orders  of  approximation  cannot  be  answered  exactly  because  of  the  lack 
of  analytic  solutions  with  which  to  compare.  The  effectiveness  of  the  approach 
will  be  established  by  frequent  use,  which  will  gradually  produce  a  qualitative 
feel  for  its  accuracy.  A  small  start  has  been  made  in  this  direction.  Several 
cases  have  been  run  to  illustrate  various  features  of  the  method  and  to  compare 
it  with  the  first  order  method  of  reference  1. 

The  calculated  examples  are  of  two  types.  The  first  compares  various 
options  of  the  program  to  determine  their  relative  importance  and  to  establish 
tentatively  the  proper  use  of  the  method.  The  second- type  illustrates  the 
scope  of  the  method  by  presenting  results  for  realistic  examples. 

A  series  of  cases  was  run  with  the  edge  vortex  formulas  of  section  6.0 
and  the  global  vorticity  feature  of  section  7.2  to  establish  the  relative 
importance  of  spanwise  and  streamwise  edge  vortices  on  a  panel.  The  body  used 
for  this  study  is  the  swept  wing  of  Figure  2a.  Calculations  were  performed  with 
the  first-order  method  of  reference  1  and  with  the  present  method  suppressing 
both  curvature  and  source-derivative  effects.  Differences  between  the  two  cases 
are  due  solely  to  edge-vortex  effects.  If  the  edge-vortex  terms  of  section  6.0 
are  used  along  all  N-lines  the  present  method  and  the  first-order  method  have 
the  same  trailing  vorticies.  If  either  the  piecewise-constant  spanwise  vort’city 
option  is  used  or  the  piecewise-linear  spanwise  vorticity  option  is  used  with 
the  global  vorticity  feature  of  section  7.2,  then  neither  method  has  vortices 
along  spanwise  panel  edges.  If  the  features  of  both  the  last  two  sentences  are 
used,  then  the  two  programs  should  give  identical  answers.  This  was  verified 
explicitly.  If  the  N-line  edge  vortices  are  retained  but  the  piecewise-linear 
spanwise  vorticity  option  is  used  without  the  feature  of  section  7.2,  the 
first-order  method  has  vortices  along  spanwise  panel  edges  while  the  present 
method  does  not.  A  comparison  of  the  calculated  results  then  indicates  the 
importance  of  the  spanwise  edge  vortices.  The  two  span-load  distributions  are 
essentially  Identical  as  are  the  pressure  distributions  except  for  the  peak 
velocities,  which  are  slightly  different.  It  is  concluded  that  the  use  of  the 
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feature  of  section  7.2  is  desirable  because  it  offers  some  improvement,  but 
that  the  level  of  discrepancy  does  not  warrant  special  edge-vortex  formulas 
for  spanwise  panel  edges  for  curved-panel  cases  where  exact  cancellation  of 
these  vortices  cannot  be  attained. 

The  study  of  the  previous  paragraph  was  continued  by  using  the  present 
method  with  an  edge  vortex  on  the  tip  N-line  only  and  comparing  its  results  to 
those  obtained  using  edge-vortices  on  all  N-lines.  To  minimize  the  strengths  of 
the  interior  edge  vortices,  the  pi ecewise-1 inear  spanwise  vorticity  option  was 
used.  Figure  3  compares  the  two  spanwise  distributions  of  section  lift  coef¬ 
ficients.  It  can  be  seen  that  the  neglect  of  trailing  vortices  on  the  interior 
N-lines  yields  a  somewhat  wiggly  lift  distribution  of  reduced  level.  Accordingly, 
it  was  decided  to  include  these  vortices  in  all  future  calculations. 

As  described  in  section  8.0,  the  Kutta  condition  in  the  first-order  method 
has  been  applied  by  requiring  equal  pressures  at  the  control  points  nearest  the 
trailing  edge  on  the  upper  and  lower  surfaces,  i.e.,  the  equal-pressure  condi¬ 
tion  is  applied  a  finite  distance  forward  of  the  trailing  edge.  On  supercritical 
wings,  for  which  pressure  gradients  are  large  in  the  trailing  edge  regions,  this 
can  lead  to  an  underprediction  of  lift.  As  a  remedy,  a  new  form  of  the  Kutta 
condition  has  been  developed  in  which  upper-  and  lower-surface  pressures  are 
separately  extrapolated  to  the  trailing  edge,  and  these  extrapolated  values  are 
set  equal.  Sample  results  are  shown  in  Figure  4  for  the  supercritical  wing  whose 
planform  is  shown  in  Figure  2b  and  whose  airfoil  section  is  as  given  in  Figure  4. 
As  is  evident  in  Figure  4a,  the  requirement  of  pressure  equality  forward  of  the 
trailing  edge  pulls  the  pressure  curve  inside  the  more  correct  one.  The  result¬ 
ing  loss  of  lift  is  more  easily  seen  in  the  spanwise  distribution  of  Figure  4b 
and  in  the  values  of  total  lift.  Thus  the  extrapolated  Kutta  condition  is  used 
exclusively  in  the  present  method. 

To  compare  the  present  method  with  the  first-order  method  of  reference  1, 
calculations  were  performed  for  the  conventional  wing  whose  planform  is  shown 
in  Figure  2a  and  whose  airfoil  section  is  given  in  Figure  5  and  also  for  the 
wing-fuselage  shown  in  Figure  6.  Results  for  the  clean  wing  are  shown  in  Fig¬ 
ure  5.  As  expected  from  two-dimensional  experience,  there  is  not  a  great  deal 
of  difference  between  the  two  methods  for  such  a  simple  geometry.  Nevertheless, 
some  change  in  pressure  distribution  Is  evident  In  Figure  5a.  The  most  sig¬ 
nificant  difference  is  in  the  spanwise  lift  distribution  shown  in  Figure  5b 
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where  It  is  seen  that  the  higher-order  method  gives  higher  tip  loadings  and 
lower  root  loadings. 

Despite  its  greater  geometrical  complexity,  the  results  for  the  wing-fuselage 
are  not  very  different  for  the  two  programs,  as  shown  in  Figure  7.  The  higher- 
order  method  yields  a  smoother  spanwise  variation  of  section  lift  coefficient. 


(x,y,z) 


-n 


Figure  1.  A  general  curved  surface  panel  and  its  projection  in  the  tangent 
plane. 
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Figure  2.  Planforms  and  panel  distributions  for  two  wings 
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options  for  the  trailing  vorticity. 


(b)  Spanwise  variations  of  section  lift  coefficient. 


Figure  4 


Calculated  results  for  the  supercritical  wing  with  and  without  the 
extrapolated  Kutta  condition. 
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(a)  Chordwlse  pressure  distributions  near  the  root 


•»  ■ 

(b)  Spanwlse  variation  of  section  lift  coefficient 
Figure  7.  Calculated  results  for  a  wing-fuselage  combination 
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APPENDIX  A 

A  SURFACE  GEOMETRY  FITTING  PROCEDURE  BASED  ON  BICUBIC  SPLINES 

The  original  -le  st-square  fitting  procedure  for  P,  Q  and  R,  as  described 
in  section  7.2  of  reference  2,  proved  satisfactory  for  many  geometries.  However, 
it  broke  down  for  large  panel  aspect  ratios,  as  for  example  occur  near  wing 
leading  edges.  Several  modifications  were  attempted,  including  weighted  least 
squares  and  double-precision  arithmetic,  but  none  proved  satisfactory  for  every 
application.  Accordingly,  it  was  decided  to  try  a  completely  new  geometric 
procedure  to  replace  that  of  section  7.0  of  reference  2.  This  approach  effec¬ 
tively  decouples  the  geometric  calculation  from  the  fluid  dynamic  one,  and  it  has 
proved  more  accurate  and  versatile  than  any  of  its  predecessors. 

Very  elaborate  geometry  fitting  procedures  based  on  parametric  bicubic 
splines  have  been  developed  at  Douglas  Aircraft  Company  over  many  years.  A 
description  of  this  technique  is  beyond  the  scope  of  the  present  report.  A 
survey  is  contained  in  reference  8.  In  the  present  application  the  method  is 
considered  a  "black  box,"  although  several  minor  changes  had  to  be  made. 

The  points  defining  the  body  are  input  in  the  usual  way.  Each  panel  is 
fitted  by  a  bicubic  surface  in  terms  of  two  parameters  u  and  v  that  vary 
from  0  to  1  over  the  panel.  (The  panel  is  the  unit  square  in  parameter  space.) 
This  permits  the  well-known  procedures  of  reference  9  to  be  used  as  follows. 

Let  a  point  (x,y,z)  of  the  panel  be  represented  as  a  vector 

x  =  xT  +  yj  +  z£ 

The  parametric  cubic  fit  then  yields 

x  =  x(u,v) 

These  expressions  may  be  differentiated  analytically  to  give 


as  functions  of  u  and  v.  The  vectors  xu  and  xy  are  tangent  to  the 
curves  v  =  constant  and  u  =  constant,  respectively,  and  thus  lie  in  the 
surface  although  they  are  not  perpendicular. 


(A.l) 
(A. 2) 
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The  point  corresponding  to  u  =  v  =  H  is  in  the  "center"  of  the  panel  in 
some  sense.  It  is  selected  as  the  control  point  and  origin  of  coordinates  of  the 
flat  projected  panel.  The  derivatives  of  (A. 3)  are  evaluated  there,  and  in  all 
that  follows  x  and  its  derivatives  are  assumed  to  be  those  at  u  =  v  = 


The  unit  normal  vector  to  the  panel,  which  is  also  the  unit  vector  along 
the  c  axis  of  panel  coordinates  is 


-► 

n  = 


(A. 4) 


where  the  sign  is  selected  to  give  an  outward  normal.  The  unit  vector  along 
the  s  axis  of  panel  coordinates  is  taken  tangent  to  the  v  =  constant  curve 
which  nearly  parallels  the  N-lines, 


(A. 5) 


Thus  the  unit  vector  along  the  n  axis  of  panel  coordinates  is 


(A. 6) 


The  components  of  the  three  unit  vectors  thus  obtained  are  the  transformation 
matrix.  Compare  equation  (7.2.10)  of  reference  1. 

Now  define 

h  *  u  -  H,  k 

and  consider  the  Maclaurin  series  for  £,  n. 

They  have  the  form 

5  =  Ah  +  Bk  +  (second  order) 
n  =  Ch  +  Dk  +  (second  order) 
c  =  *s(eh2  +  2fhk  +  gk2)  +  (third  order) 

There  are  no  constant  terms  in  (A. 8),  because  the  origin  of  panel  coordinates 
corresponds  to  h  =  k  =  0.  Furthermore,  since  the  £n  plane  Is  tanyont  to 
the  surface  at  the  origin  (£e  is  the  normal  vector),  the  series  for  c  has 
no  linear  terms.  Reference  9  gives  the  coefficients  of  equations  (A. 8)  as 


=  v  -H  (A. 7) 

and  e  in  terms  of  h  and  k. 

(A. 8) 
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VTe 

c  ■  K  •  I.  ■ 0 


B  =  x. 


D  =  xy  •  Je 


(A. 9) 


e  =  xuu  '  n 


f  =  V  •  * 


9  »  *w  *  " 


(A. 10) 


The  first  two  of  equations  (A. 8)  may  be  inverted  to  give 


h  =  a?  +  bn  +  (second  order) 
k  =  cc  +  dn  +  (second  order) 


(A. 11) 


where 


a=  a  • 


c  =  I  ’ 


A  =  AD  -  BC 


(A. 12) 


Equation  (A. 12)  may  be  inserted  into  the  third  equation  of  (A. 8)  to  give  the 
desired  form 


c  -  Pc2  +  2Qen  +  Rn2 


(A. 13) 


The  result  is 


P  =  *s[ea2  +  2fac  +  gc2] 

Q  =  *s[eab  +  f(ad  +  be)  +  ged] 
R  =  4(eb  +  2fbd  +  gd2] 


(A. 14) 


For  generality  c  has  been  included  in  (A. 14),  but  in  the  present  application 
it  is  zero,  which  simplifies  (A. 14). 

It  remains  to  compute  corner  points  in  panel  coordinates.  The  four  input 
points  bounding  the  panel  are  transformed  into  panel  coordinates  to  obtain 
(C|,  c£)  k  =  1,  2,  3,  4.  They  are  projected  into  the  plane  by  simply 

ignoring  c*.  Next  the  side  between  points  1  and  2  is  rotated  to  make  =  n2- 
The  midpoint  and  length  of  the  side  are,  respectively. 
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e-JsUf  +  ef).  n  =  >s(nf  +  n|) 

d  =  ^(5f  -  qf  +  (nf  -  T,|)2 

Then  the  final  corner  point  coordinates  are 

n-,  =  n2  =  n 


(A. 15) 


(A. 16) 


A  similar  calculation  is  performed  for  the  side  between  the  points  3  and  4. 

It  should  be  noted  that  the  underlying  parametric  cubic  geometry  routine 
uses  the  surrounding  input  points  to  generate  the  fit  to  a  panel.  The  routine 
considers  only  points  on  the  same  section,  and  thus  slightly  different  results 
can  be  obtained  depending  on  how  the  body  is  sectioned.  For  fitting  purposes 
the  wake  is  considered  a  separate  section,  so  that  the  routine  does  not  try  to 
fit  around  the  trailing  edge.  On  the  semi -infinite  last-wake  panel  the  deriv¬ 
atives  P  and  Q  are  set  equal  to  zero,  so  that  the  panel  has  straight  gener¬ 
ators  in  the  stream  direction,  but  R,  the  spanwise  second  derivative,  may  be 
nonzero. 


The  parametric  cubic  procedure  gives  smooth  fits  even  in  extreme  cases. 

By  way  of  illustration  the  airfoil  sections  on  a  wing  were  defined  by  only  four 
points,  leading  and  trailing  edges  and  upper  and  lower  surface  maximum-thickness 
points.  Thickness  was  20%  at  the  root  and  mid-semispan  and  10%  at  the  tip. 
Results  are  shown  in  Figure  A.l.  Remarkably,  the  "diamond"  section  shapes  have 
been  fiited  with  reasonable,  smooth  curves,  as  opposed,  for  example,  to  what  a 
parabolic  fit  would  have  done.  It  should  be  noticed  that  the  leading  edge  is 
rounded  and  the  trailing  edge,  from  which  the  wake  issues,  is  sharp  as  required. 
The  fact  that  the  trai ling-edge  angles  are  somewhat  large  is  due  to  the  absurdly 
inadequate  point  number. 
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Parametric  cubic  representation  of  a  wing  with  four  panels  around 
each  section. 


Unclassified 


SECURITY  CLASSIFICATION  OF  THIS  PAGE  (When  Pete  Entered) _ 

REPORT  DOCUMENTATION  PAGE  befoI4DcSSJletog  form 

1.  REPORT  NUMBER  ~~  ]2.  GOVT  ACCESSION  NO.  3-  RECIPIENT'S  CATALOG  NUMBER 


1.  REPORT  NUMBER 

NADC- 79277-60 

4.  TITLE  (end  Subtitle) 


-(Uzd 


An  Improved  Higher-Order  Panel  Method  for 
Three-Dimensional  Lifting  Potential  Flow 


17.  author^; 


5.  TYPE  OF  REPORT  «  PERIOD  COVERED 


Final  Report 


PERFORMING  ORG.  REPORT  NUMBI 

MDC  J2162 

I.  CONTRACT  OR  GRANT  *UMSERf,J 


John  L.  Hess 


N62269-80-C-0228 


9.  PERFORMING  ORGANIZATION  NAME  AND  ADDRESS 

Douglas  Aircraft  Company 
3855  Lakewood  Boulevard 
Long  Beach,  California  90846 

11.  CONTROLLING  OFFICE  NAME  AND  ADDRESS 

Naval  Air  Development  Center 
Warminster,  PA  18904  (6052) 


to.  PROGRAM  ELEMENT.  PROJECT.  TASK 
AREA  *  WORK  UNIT  NUMBERS 


N62269-80-65-00495 


12.  REPORT  DATE 

December  1981 

13.  NUMBER  OF  PAGES 


14.  MONITORING  AGENCY  NAME  A  ADDRESS  (U  different  from  Controlling  Office)  13.  SECURITY  CLASS,  (ot  (Ms  report) 

Unclassified 

'  lift.  OECLASSIFICATION/ DOWN  GRADING 
SCHEDULE 

1«.  DISTRIBUTION  STATEMENT  (of  thle  Report) 

Approved  for  Public  Release;  Distribution  Unlimited 


[17.  DISTRIBUTION  STATEMENT  (of  the  ebetrect  entered  In  Block  20,  It  different  from  Report) 


ie.  supplementary  notes 


1 19.  KEY  WOROS  (Continue  on  rereree  elde  It  neceeeery  end  Identity  by  block  number) 


Aerodynamics 
Computer  Program 
Flow  Field 
Fluid  Dynamics 


Higher  Order 
Numerical  Analysis 
Panel  Method 
Potential  Flow 


Source  Density 
Surface  Singularity 
Three-Dimensional  Flow 


.  ABSTRACT  (Continue  on  reeeree  tide  It  neceeeery  md  Identity  by  block  number) 

kA  previous  report  described  a  higher-order  three-dimensional  panel  method 
that  represents  the  body  about  which  flow  Is  to  be  computed  by  means  of 
curved  four-sided  surface  panels  having  linearly  varying  source  and  vortlclty 
distributions.  The  method  of  accounting  for  lift  was  Incomplete,  and  the 
main  purpose  of  the  present  work  was  to  remedy  this  defect.  A  number  of 
other  modifications  to  the  method  were  also  made  to  Improve  the  efficiency 
and  accuracy  of  the  method. 


FORM  im 
I  JAN  7*  WJ 


EDITION  OF  I  NOV  •>  I*  OBSOLETE 
S/N  0102*014*  SS0I  I 


_ Unclassified _ 

SECURITY  CLASSIFICATION  OF  THIS  PAGE  I 


CIJMITY  CLASSIFICATION  of  this  PAOEl'H'h«n  D«t« 


-The  modifications  documented  here  are:  formulas  for  the  edge-yortex  influ¬ 
ences,  which  previously  had  been  neglected;  new  surface  vorti city  formulas 
that  express  its  influences  in  terms  of  source  Influences;  a  modified  global 
vorti city  algorithm  to  Improve  continuity  over  the  surface,  and;  an 
extrapolated  Kutta  condition.,.. 


_ Unclassified _ 

SSCUHITT  CLASSIFICATION  OF  THIS  NAOer*fc*n  Dmtm  Knt*n4) 


DISTRIBUTION  LIST 


Report  No.  NADC-79277-60 
AIRTASK  No.  A320320D/001A/1R023-02-000 

No.  of 
Copies 

1UWAH  (AIR-95QD) . * 

(2  for  retention) 

(1  for  AIR-320D) 

(1  for  AIR-5301) 

HAWFMCEH,  Chino  Lake,  CA . 1 

NAVAIRPROPCEH,  Trenton,  RJ.  .  . . 1 

DTNSKDC,  Bethesda,  MD  (Attn:  Dr.  H.  Chaplin)  . . 1 

ONR,  Arlington,  VA  (Attn:  R.  Whitehead) . 1 

HAVPGSCOL,  Monterey,  CA  (Attn:  M.  Platzer)  ...............  1 

NASA,  Aims  Research  Center,  Moffett  Field,  CA. .............  2 

(1  for  D.  Rickey) 

(1  for  V.  Deckert) 

RASA,  Langley  Research  Center,  Hampton,  VA  (Attn:  R.  Kargason) . 1 

RASA,  Levis  Research  Center,  Cleveland,  (9 . 1 

Wright-Pat terson  AFB,  Dayton,  OR.  ....................  2 

Cl  for  Flight  Dynamics  Lab) 

(1  for  Aeronautical  Systems  Division) 

The  Pentagon,  Washington,  DC  (Attn:  R.  Sievert) . 1 

fJ.S.  Army  Aviation  Systems  Command,  St.  Louis,  MD . 1 

|[.S.  Army  Research  Office,  Durham,  MC . 1 

l»TIC,  Alexandria,  VA . 12 

Boeing  Company,  Seattle,  WA  (Attn:  E.  Omar) . 1 

LTV  Aerospace  Corporation,  Dallas,  TZ...  ...............  2 

(1  for  T.  Beatty) 

(1  for  W.  Slnpkln) 

Rockwell  International,  Columbus,  OH  (Attn:  9.  Palmer)  .........  1 

General  Dynamics  Corporation,  Ft.  Worth,  TZ  (Attn:  V.  Folley) .  .....  1 

Nielson  Engineering,  Mountain  View,  CA  (Attn:  S.  Spangler) . 1 

Onlv.  of  Tennessee,  Space  Inst.,  Tullahoma,  TN  (Attn:  V.  Jacobs)  ....  1 

Lockheed-Callfomla  Co.,  Burbank,  CA  (Attn:  T.  Chin)  . . 1 

Northrop  Corporation,  Hawthorne,  CA  (Attn:  P.  Wooler).  .........  1 

Grumman  Aerospace  Corp. ,  Bethpage,  LX,  HT  (Attn:  S.  Kalamaris).  .....  1 

Royal  Aeronautical  Establishment,  Bedford,  England  (Attn:  A.  Woodfleld).  1 
Fair child- Republic,  Corporation,  Farmlngdale,  LX,  ITT.  ..........  1 

Calspan,  Buffalo,  NT.  ..........................  1 

McDonnell  Douglas  Corp.,  St.  Louis,  MD  (Attn:  Dr.  D.  Ko tans  ley)  .....  1 

f/STOL  Consultant,  Newport  Mews,  VA  (Attn:  R.  Rida).  ..........  1 

Georgia  Inst.,  of  Technology,  Atlanta,  GA  (Attn:  Dr.  I.  McMahon)  ....  1 

Penn  Stats  Unlv.,  Dnlv.  Park,  PA  (Attn:  Prof.  B.  W.  McCormick)  .....  1 


